library(ggplot2)
library(tidytext)
library(tidyverse)
library(ggplot2)
library(magrittr)
library(data.table)
library(dplyr)
library(foreign)
library(haven)
library(estimatr)
library(stargazer)
library(quanteda)
library(ri2)
library(gtable)
library(gridExtra)
library(dotwhisker)
library(gapminder)
library(broom)
library(readxl)
library(ggalluvial)
library(ggeasy)
library(modelsummary)
library(sandwich)
library(plm)
library(lmtest)

rm(list = ls())

setwd(dirname(rstudioapi::getSourceEditorContext()$path)) #sets WD to the location of this code file 

data <- read.csv("cleaned_data.csv")


#Filter treatment to variables that do not reveal identity of respondents:



# Start of Main Tables  ---------------------------------------------------

# Table 1 -----------------------------------------------------------------

trumpvote <- lm(trump_vote ~  treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
trumpvoteint <- lm(trump_vote ~  treatment + party + College + trumpfavor + treatment:trumpfavor + justice_partisan + Male + Age, data = data)

se1 <- sqrt(diag(vcovHC(trumpvote, type = "HC1")))
se2 <- sqrt(diag(vcovHC(trumpvoteint, type = "HC1")))

stargazer(trumpvote, trumpvoteint, type = "latex", label = "table1", 
          title = "Legal Rhetoric Reduces Intention of Voting for Trump, But Only among Non-Trump Supporters",
          star.cutoffs = c(.05, .01, .001),
          keep = c("treatmenttrump", "treatmentsmith", "trumpfavor"),
          dep.var.labels = c("Trump Primary Vote"),
          covariate.labels = c("Partisan Rhetoric", "Legal Rhetoric", "Trump Favorability", "Partisan Rhetoric:Trump Favorability", "Legal Rhetoric:Trump Favorability"),
          se = list(se1, se2),
          keep.stat = c("n", "adj.rsq"))

trumpvoteint <- feols(trump_vote ~  treatment + party + College + justice_partisan + Male + Age,
                      data = data, split =~ trumpfavor)

# Figure 1 -----------------------------------------------------------------

#trump_post_n is Trump.Approval on a 7 point scale negative to positive scale 
#trump_pre_n is Trump.Favorability on very unfavorable to very favorable scale 

coefs <- c(
  "treatmentsmith" = "Treatment Smith",
  "treatmenttrump" = "Treatment Trump",
  "treatmentsmith:trumpfavor" = "Smith × Trump Favor",
  "treatmenttrump:trumpfavor" = "Trump × Trump Favor")

trump_levels <- c("Very unfavorable","Somewhat unfavorable", "Neither favorable nor unfavorable",
                  "I don't know this person.", "Somewhat favorable", "Very favorable")
data$Trump.Favorability <- factor(data$Trump.Favorability, levels = trump_levels)

temp <- filter(data, Trump.Favorability != "I don't know this person.")

trumpvote5 <- feols(trump_vote ~ treatment + party + College + treatment + justice_partisan + Male + Age, data = temp, subset =~ Trump.Favorability == "Very unfavorable", vcov = "HC1")
trumpvote4 <- feols(trump_vote ~  treatment + party + College + treatment + justice_partisan + Male + Age, data = temp, subset =~ Trump.Favorability == "Somewhat unfavorable", vcov = "HC1")
trumpvote3 <- feols(trump_vote ~  treatment + party + College + treatment + justice_partisan + Male + Age, data = temp, subset =~ Trump.Favorability == "Neither favorable nor unfavorable", vcov = "HC1")
trumpvote2 <- feols(trump_vote ~  treatment + party + College + treatment + justice_partisan + Male + Age, data = temp, subset =~ Trump.Favorability == "Somewhat favorable", vcov = "HC1")
trumpvote1 <- feols(trump_vote ~  treatment + party + College + treatment + justice_partisan + Male + Age, data = temp, subset =~ Trump.Favorability == "Very favorable", vcov = "HC1")

models <- list(
  "Very favorable" = trumpvote1,
  "Somewhat favorable" = trumpvote2,
  "Neither favorable nor unfavorable" = trumpvote3, 
  "Somewhat unfavorable" = trumpvote4,
  "Very unfavorable" = trumpvote5
)

coef_plot1 <- modelplot(models, coef_map = coefs) + geom_vline(xintercept = 0, linetype = "dashed") + 
  theme_bw(base_size = 18) + 
  scale_y_discrete(labels = c("Legal Rhetoric", "Partisan Rhetoric"))  + 
  scale_color_discrete(c(" Pre-Treatment \n Trump Attitudes"))  + guides(color = guide_legend(reverse = TRUE))

coef_plot1

ggsave("coef_plot1.pdf")

# Table 2 -----------------------------------------------------------------

# Components of Prosecution Index, Full Sample  -----------------------------------------
facts <- lm(facts_numeric ~  treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
justified <- lm(justified_numeric ~  treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
norm <- lm(norm_numeric ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
pros_index <- lm(pros_index ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)

se1 <- sqrt(diag(vcovHC(facts, type = "HC1")))
se2 <- sqrt(diag(vcovHC(justified, type = "HC1")))
se3 <- sqrt(diag(vcovHC(norm, type = "HC1")))
se4 <- sqrt(diag(vcovHC(pros_index, type = "HC1")))

stargazer(facts, justified, norm, pros_index, type = "latex", label = "table2", 
          title = "Legal Rhetoric Increases Perceptions of the Prosecution as Apolitical and Consistent with Democratic Norms",
          star.cutoffs = c(.05, .01, .001),
          keep = c("treatmenttrump", "treatmentsmith"),
          dep.var.labels = c("Factually True", "Not Politically Motivated", "Consistent With Dem Norms", "Index of Support"),
          covariate.labels = c("Partisan Rhetoric", "Legal Rhetoric"),
          se = list(se1, se2, se3, se4),
          keep.stat = c("n", "adj.rsq"))

# Table 3 -----------------------------------------------------------------

pros_index <- lm(smith_change ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
pros_CATE <- lm(smith_change ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age + treatment:trumpfavor, data = data)

se1 <- sqrt(diag(vcovHC(pros_index, type = "HC2")))
se2 <- sqrt(diag(vcovHC(pros_CATE, type = "HC2")))

stargazer(pros_index, pros_CATE, type = "latex", label = "table3", 
          title = "Legal Rhetoric Strongly Reduces Support for the Prosecutor among Trump Supporters",
          star.cutoffs = c(.05, .01, .001),
          keep = c("treatmenttrump", "treatmentsmith"),
          dep.var.labels = c("Approval of Prosecutor", "Approval of Prosecutor"),
          covariate.labels = c("Partisan Rhetoric", "Legal Rhetoric", "Partisan Rhetoric*Trump Favorability", "Legal Rhetoric:Trump Favorability"),
          se = list(se1, se2),
          keep.stat = c("n", "adj.rsq"))

# Table 5  ----------------------------------------------------------------

norm1 <- lm(direct_index~  treatment + treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data)
norm2 <- lm(indirect_index ~  treatment + trumpfavor + party + College + justice_partisan + Male + Age, data = data)
norm3 <- lm(affpol ~  treatment + trumpfavor + party + College + justice_partisan + Male + Age, data = data)

se1 <- sqrt(diag(vcovHC(norm1, type = "HC1")))
se2 <- sqrt(diag(vcovHC(norm2, type = "HC1")))
se3 <- sqrt(diag(vcovHC(norm3, type = "HC1")))

stargazer(norm1, norm2, norm3, type = "latex", label = "table5", 
          title = "Partisan and Legal Rhetoric Have No Overall Effect on Support for Democratic Norms or Affective Polarization",
          star.cutoffs = c(.05, .01, .001),
          keep = c("treatmentsmith", "treatmenttrump"),
          dep.var.labels = c("Support for Direct Retaliation", "Support for Indirect Retaliation", "Affective Polarization"),
          covariate.labels = c("Partisan Rhetoric", "Legal Rhetoric"),
          se = list(se1, se2, se3),
          keep.stat = c("n", "adj.rsq"))


norm1I <- feols(direct_index ~ treatment + College + trumpfavor + justice_partisan + Male + Age, data = data, subset =~ party == "independent")
norm1L <- feols(direct_index ~ treatment + College + trumpfavor + justice_partisan + Male + Age, data = data, subset =~ party == "lean")
norm1R <- feols(direct_index ~ treatment + College + trumpfavor + justice_partisan + Male + Age, data = data, subset =~ party == "republican")

coef_plot2 <- modelplot(list("Strong Rep" = norm1R, "Lean Rep" = norm1L, "Independents" = norm1I), coef_map = coefs) + geom_vline(xintercept = 0, linetype = "dashed") + 
  theme_bw(base_size = 18) + 
  scale_y_discrete(labels = c("Legal Rhetoric", "Partisan Rhetoric"))  + 
  scale_color_discrete(c("Party ID"))  + guides(color = guide_legend(reverse = TRUE))

coef_plot2

ggsave("coef_plot2.pdf")

# APPENDIX STARTS HERE  ---------------------------------------------------

#Descriptive Statistics: Note many of these are re-labelled in the text 
stargazer(data, summary = TRUE)

#Knowledge of democratic norms
prop.table(table(data$Democracy.1))
prop.table(table(data$Democracy.2))


#Create CNBC perception plot

data$CNBC.Perception <- as.factor(data$CNBC.Perception)
data$CNBC.Perception <- factor(data$CNBC.Perception, levels = c("Very Democrat", "Somewhat Democrat", "Neither Republican nor Democrat",
                                                                "Somewhat Republican", "Very Republican"))

freq_table <- as.data.frame(table(data$CNBC.Perception))
colnames(freq_table) <- c("category", "count")
freq_table$percentage <- (freq_table$count / sum(freq_table$count)) * 100

ggplot(freq_table, aes(x = category, y = percentage)) +
  geom_bar(stat = "identity") +
  ylim(0, 50) +  # Set y-axis limits to 0-100 for percentage
  labs(title = "", 
       x = "Category", 
       y = "Percentage (%)") +
  theme_minimal()

#Create Justice Department perception Plot

data$Non_partisan <- factor(data$Non_partisan, levels = c("Not at all partisan", 
                                                          "Not very partisan", "Somewhat partisan", "Very partisan", "Extremely partisan"))

freq_table <- as.data.frame(table(data$Non_partisan))
colnames(freq_table) <- c("category", "count")
freq_table$percentage <- (freq_table$count / sum(freq_table$count)) * 100

ggplot(freq_table, aes(x = category, y = percentage)) +
  geom_bar(stat = "identity") +
  ylim(0, 50) +  # Set y-axis limits to 0-100 for percentage
  labs(title = "", 
       x = "Category", 
       y = "Percentage (%)") +
  theme_minimal()

ggsave("partisan.pdf")


# Create Pre-Treatment Trump Favor Plot ----------------------------------


freq_table <- as.data.frame(table(data$Trump.Favorability))
colnames(freq_table) <- c("category", "count")
freq_table$percentage <- (freq_table$count / sum(freq_table$count)) * 100

ggplot(freq_table, aes(x = category, y = percentage)) +
  geom_bar(stat = "identity") +
  ylim(0, 50) +  # Set y-axis limits to 0-100 for percentage
  labs(title = "", 
       x = "Category", 
       y = "Percentage (%)") +
  theme_minimal()


ggsave("trump.pdf")


#Explore treatment-specific attention checks 

data %<>% filter(!is.na(Placebo.Attn.Check) | !is.na(nonpartisan_check) | !is.na(partisan_check))

prop.table(table(data$Placebo.Attn.Check, data$treatment), 2)
prop.table(table(data$nonpartisan_check, data$treatment), 2)
prop.table(table(data$partisan_check, data$treatment), 2)


# Multiple-Comparison Corrections, Main Outcomes  ----------------------------------------

# Benjamini-Hochberg Correction -------------------------------------------

#We consider our main 6 outcome variables, running in lm_robust to make extracting P values easier 

pol1 <- lm_robust(pros_index ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data, se_type = "HC2")
pol2 <- lm_robust(smith_change ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data, se_type = "HC2")
pol3 <- lm_robust(trump_vote ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data, se_type = "HC2")
norm1 <- lm_robust(direct_index ~ treatment + party + trumpfavor + treatment + justice_partisan + Male + Age, data = data, se_type = "HC2")
norm2 <- lm_robust(indirect_index ~ treatment + party + trumpfavor + treatment + justice_partisan + Male + Age, data = data, se_type = "HC2")
norm3 <- lm_robust(affpol ~ treatment + party + College + trumpfavor + justice_partisan + Male + Age, data = data, se_type = "HC2")

pvalues1 <- c(pol1$p.value[2], pol2$p.value[2], pol3$p.value[2], norm1$p.value[2], norm2$p.value[2], norm3$p.value[2])
pvalues2 <- c(pol1$p.value[3], pol2$p.value[3], pol3$p.value[3], norm1$p.value[3], norm2$p.value[3], norm3$p.value[3])
format(pvalues1, scientific = F, digits = 3)
format(pvalues2, scientific = F, digits = 1)

ranks1 <-rank(pvalues1, ties.method = "last")
ranks2 <-rank(pvalues2, ties.method = "last")

fdrs1 <-p.adjust(pvalues1, method="BH")
fdrs2 <-p.adjust(pvalues2, method="BH")

format(fdrs1, scientific = F, digits = 3)
format(fdrs2, scientific = F, digits = 3)





